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Abstract 

We formulate a fully self-consistent, microscopic model to study the retardation 
and correlation effects of the barrier within a Josephson junction. The junction 
is described by a series of planes, with electronic correlation included through a 
local self energy for each plane. We calculate current-phase relationships for var- 
ious junctions, which include non-magnetic impurities in the barrier region, or an 
interfacial scattering potential. Our results indicate that the linear response of the 
supercurrent to phase across the barrier region is a good, but not exact indicator of 
the critical current. Our calculations of the local density of states show the current- 
carrying Andreev bound states and their energy evolution with the phase difference 
across the junction. 

We calculate the figure of merit for a Josephson junction, which is the product 
of the critical current, /cj arid the normal state resistance, i?7v, for junctions with 
different barrier materials. The normal state resistance is calculated using the Kubo 
formula, for a system with zero current flow and no superconducting order. Semi- 
classical calculations would predict that these two quantities are determined by the 
transmission probabilities of electrons in such a way that the product is constant 
for a given superconductor at fixed temperature. Our self-consistent solutions for 
different types of barrier indicate that this is not the case. We suggest some forms 
of barrier which could increase the IcRn product, and hence improve the frequency 
response of a Josephson device. 

Typeset using REVI^ 



I. INTRODUCTION 



Some of the most successful and most promising electronics applications of superconductors 
involve Josephson junctions, where two superconducting regions are coupled through a barrier 
region, made from a non-superconducting material. The drive to make electronic devices using 
Josephson junctions has been motivated by their naturally high operating frequencies, far in excess 
of the speeds obtainable in standard silicon technology. The operational frequency of an ideal 
Josephson junction using rapid single flux quantum (RSFQ) logic is /c-Ratc/^, where Ic is the critical 
current, and Rn is the normal state resistance of the junctionl^hl. Hence a maximum value of the 
product, IcRn leads to optimal performance. At present, low temperature superconductors have 
been used to produce junctions with IcRn products up to ImV, and operational speeds reaching 770 
ghJ'I, while junctions of high-temperature superconductors have achieved IcRn products reaching 
from 1 mV to ~ 20 mvli^. The choice of barrier material strongly affects both Ic and Rn (typically 
poor conductors increase Rn but reduce Ic) so it is appropriate to study in detail what properties 
of the junction modify /c, Rn, and IcRn- Indeed, it has been suggested^, that the maximum value 
of the product would be reached when the barrier material is close to a metal-insulator transition, 
so that the system is near the cross-over from an SNS to an SIS junction. In order to investigate 
such non-trivial barrier materials, we develop a microscopic model of a Josephson junction, which 
self-consistently incorporates the dynamical correlation effects of the electrons in calculating the 
conductance of the junction and the critical current. 

Our model is also appropriate for ballistic junctions that have relatively pure barrier regions so 
that the extent of the barrier is smaller than its mean-free-path but larger than its proximity-effect- 
induced superconducting correlation length. Much work has progressed on these systems recently, 
with a concentration on Niobium-Indium Arsenide-Niobium junctions or Niobium- Silicon- Niobium 
junctions. Recent work includes an examination of subgap structures and current deficitJiil, quasi- 
particle reflection effects and spikes in the conductanceEl, effects of multiple Andreev reflection and 
noise in stacked junctionslll, and an investigation of a tunable junction that can be altered from an 
ordinary junction to a 7r-junction by driving the barrier into a nonequilibrium state via a transverse 
electrical curren t0. The above-mentioned experimental and theoretical work has concentrated on 
nonequilibrium effects that occur at finite voltages. This contribution will address only equilibrium 
and linear- response properties of the microscopic model for Josephson junctions, but our work can 
be extended to examine such nonequilibrium effects as well. 
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Josephson junctions consist of two superconducting regions coupled through an intermediate 
region which is not naturally superconducting. The main characteristics of a Josephson junction 
can be understood by just considering two superconductors coupled together with a single, energy 
independent transmissivity parameter. The dc Josephson effect, of a direct supercurrent through 
the junction at zero external voltage, and the ac Josephson effect of an oscillating, dissipative cur- 
rent at finite voltage were both predicted from such a simple mo del0. Further details of the I-V 
characteristics, such as excess current and sub-harmonic gap structure due to Andreev bound states 
were explainedEl'B, by matching boundary conditions across a step-function in the superconducting 
potential, and including an adjustable scattering potential at the interfaces0'Ei. Analytic calcula- 
tions have been performed to provide current-phase relationshipslli^, representing the barrier as a 
single scattering potential with no spatial extent. In this contribution we compare these traditional 
approaches with a more microscopic model. 

Originally the term Josephson junction referred to superconductor - insulator - superconductor 
(SIS) tunnel junctions, but here we use the term in its common usage to refer to all "weak links", 
including superconductor - normal metal - superconductor (SNS) junctions. In the latter case, the 
barrier region can occupy tens to hundreds of nanometres, as the proximity effect ensures that the 
superconductivity extends into the barrier region, decaying on the scale of the coherence length, 

= hvp/ksT when in the ballistic regime (where the coherence length is longer than the mean 
free path in the barrier). 

In this paper, we present a method for studying the effects of the barrier region on the strength 
of superconductivity it can support, and hence on the supercurrent it can maintain. Effects of 
electron correlations are incorporated within the dynamic mean field theoryiHH^, which leads to a 
local, frequency-dependent self-energy, which we allow to vary from one plane to the next, while 
assuming it to remain constant within individual planes. Hence the three-dimensional system be- 
comes inhomogeneous in the one dimension where the current flows (which we label the z-direction). 
We model the system with two sets of Nsc planes of superconducting material coupled each on one 
side to the bulk superconductor, and on the other side sandwiching Nb planes of barrier material 
as depicted in Fig. 0. The total number of planes modeled self-consistently is 2Nsc + Ni, = N. 
We note that while in this paper we concentrate on applying the method to Josephson junctions, 
our method would also be applicable to theoretical studies of correlated electrons at surfaces, single 
interfaces, or multiple interfaces, as first demonstrated by Potthoff and Noltingil. 

An important advantage of our scheme is that the different (material specific) microscopic models 



that best describe any particular material can be coupled together across the planes. Hence a 
superconductor with electron-phonon coupling described by a Holstein modelil (or even using the 
appropriate a^F in Migdal-Eliahberg theory) can be connected to a metallic region with impurities 
described by the Falicov-Kimball model (or a correlated Hubbard model), and so on. The local 
self energy for each plane is calculated independently, according to the model best-suited to the 
particular material. Once a set of local self-energies are evaluated for each plane, the Green's 
functions are calculated by finding the inverse of an infinite matrix which includes planes of bulk 
superconductor extending in the positive and negative z-direction. The inversion process, which is 
made tractable by a continued fraction representation, couples all of the different planes together, 
such that a change in the self energy on one plane affects the local interacting Green's functions on 
all other planes, particularly those nearby in real space. 

In order for a supercurrent to flow, a phase gradient must be applied to the superconductor 
and across the barrier region. The critical current, Jc, is reached, when the planes with the lowest 
superconducting order, typically at the center of the barrier region, can no longer support the 
necessary phase gradient to maintain current continuity. In our model we find the supercurrent as 
a function of phase variation across the barrier, by solving the system self-consistently at each set 
of phases. The self-consistency is crucialHl^, as the existence of a current flow affects the value of 
the superconducting order parameter, both inside and outside the barrier region. 

Section II contains a detailed description of our method, including the physical approximations 
used, and the general computer algorithm. In Section HI we analyze the current-phase relationships 
of our results in terms of the transmissivity of a barrier. The results presented in Section IV 
demonstrate the efficacy of the method to solve some simple models of Josephson junctions, using the 
Bogoliubov-de Gennes equations, where we demonstrate the effects of self-consistency, in particular 
on the superconducting order and electron density. Section V includes the results for barriers with 
impurity scattering, with a description of the calculations of normal-state resistance, and results of 
resistance and IcRn products. We conclude in Section VI with some comments on the results, and 
suggestions of other situations well-suited to our model. 



II. METHOD 
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A. Model 



Our method consists of two stages. First, we determine the properties of the bulk boundary 
regions for a uniform system. When no current flows in the bulk (homogeneous) superconductor, 
dynamical mean field theoryiH^i is employed to determine the local self energy, using the local 
approximation. When a uniform current is flowing, one must include a uniform phase variation in 
the superconducting order parameter, A{z), for the bulk system. Next, we solve the inhomogeneous 
problem self-consistently by iteration, for a number of planes, A^, coupled on either side to the 
uniform bulk solution. That is, giving the plane an index, a, and defining the central link of the 
junction to be between planes a = and a = 1, we include A*";, barrier planes, surrounded by 
Ngc self-consistently calculated superconductor planes on each side, such that N = Nf, + 2Nsc and 
planes with index, a < 1 — N/2 or a > N/2 are invariant homogeneous bulk planes (see Fig. |l|). 
Typically, we solve systems with Ngc = 30, which is significantly greater than the coherence length 
for our bulk superconductor (^ = hvp/A ^ 10a, where a is the lattice spacing). We observe that, 
except when very close to T^, the superconducting order has completely healed from its disruption 
at the interface, by the time we reach the planes at the bulk superconductor boundary (planes with 
a ^ ±N/2). 

We describe the system with the following tight-binding model: 

where c]^ and Qo- are fermionic operators which respectively create and destroy an electron of spin-cr 
in a single Wannier (tight-binding) state on the lattice site i; 

ea if i = j on plane a, 

—ta if i and j are neighboring sites on the same plane, a, 

tij 



-y/tata' if i and j are neighboring sites on consecutive planes, a and a', 
otherwise, 

with —ta the overlap or hopping integral for the a-th plane, e, the local site energy, f/j the renor- 
malized on-site, Hubbard interaction energy, Uf^ the impurity potential, and ifj = 1 if there is 
an impurity on site i, and Wi = otherwise. The magnitude of the hopping integral in the super- 
conducting region, t, is constant, and defines our energy scale for the entire system {t = 1). In 
all results that we present here, the superconducting region has an attractive Hubbard interaction. 
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Ui = —2 and no impurities {wi = for all sites on planes a < 1 ~ Nb/2 and a > Nb/2). We 
utilize the spinless Falicov-Kimball model to describe non-magnetic charge impurities within the 
barrier region. The interaction between dopant atoms and conduction electrons is attractive, and 
represented by a negative Uf^. The average impurity concentration within the barrier is given by, 
pimp = (f ) S'W^i) where A^^ (^^^ is the total number of barrier sites. We also consider systems 
where there are no impurities within the barrier, but there is a Hubbard interaction for sites within 
the barrier, Ui = Ub, that differs from Ui = —2 in the superconductor. In addition, we have included 
systems where the hopping integral, ta, differs in the barrier region, and where the interfacial planes 
{a = 1 — Nf,/2 and a = Nf,/2 only) have a non-zero local on-site potential, = Vjnt- 

B. Green's Function Calculations 

We use the matrix formulation of Nambu0 for the Green's function, G(rj, r^, iun), between two 
lattice sites and at the Matsubara frequency, iun = i'n{2n + VjksT, 

/ 

Q{Yi,Vj,iuJn) = 



G(r., rj,iujn) F{ri,rj,iuJn) 
yF{ri,rj,iuJn) -G*{ri,Vj,iuJn) 



(2) 



and the corresponding local self energy. 



E(ri,iuJn) (j)(TiJuJ„ 



g(ri, iujn) =1 I . (3) 

^*{ri,iun) -^*{ri,iun)) 

The diagonal and off-diagonal Green's functions are defined respectively as: 

G{ri,rj,iuJn) = - f drexp^iUnT) (T-rCj^{T)cl^{0)) , (4) 



F{Yi,Yj,iUn) = - (irexp(iiu;„r) (T^Cji(r)cii(O)) , (5) 
Jo 

where T,- denotes time-ordering in r and (3 = l/lksT). 

The self energies and Green's functions are coupled together through Dyson's equation, 

Q{ri,rj,iUn) = Q}^\ri,rj,iUn) + ^G}°\ri,ri,iuJn)Liri,iuJn)Q{ri,rj,iuJn), (6) 

I 

where we have included the local approximation for the self energy, S(rj, r^, = S(rj, 
[which can be relaxed if we use the Dynamical Cluster Approximation (DCA)@ii]. The non- 
interacting Green's function, Q^-^^Vi, tj, iun) is diagonal in Nambu space, with upper diagonal com- 
ponent given by: 



We emphasize that is the non-interacting Green's function and is not the effective medium of an 
equivalent atomic problem (see below for athe detailed algorithm used to solve the dynamical mean 
field theory). A major innovation in our work is to utilize an efficient hybrid real-space — momentum- 
space method for calculating the Green's functions from the set of local self-energies. We find this 
method to be much more powerful in solving systems with spatial variations or inhomogeneity, 
and it is also faster in bulk systems with current flow than a more conventional k-space integral 
technique. 

Since the stacked planes have translational symmetry within the plane, the systems that we 
study are inhomogeneous in one direction only. We choose that direction to be labeled the z-axis, 
which is also the direction of current flow through the Josephson junction. The first stage of our 
method is to convert the problem from a three-dimensional system to a one- dimensional system 
following the algorithm of Potthoff and Nolting0. We perform a Fourier transform within the 
planes to determine the planar indexed Green's functions, 

^isi^^n, ky) = ^-^^ Q{ri, Tj, iun) exp [K{xj - Xi) + ky{yj - yi)] , (8) 

Xj ,yj 

where a and (5 denote distinct planes, defined by a = Zj/a, (3 = zj/a and the summation is over all 
lattice sites, {xj,yj), within the /5-th plane. The self energy, ^^{iujn) = S(2;j,-ju;„) = S(ri,ici;„), is 
independent of the planar coordinates Xi and yi, so that Dyson's equation [Eq. (H)] becomes 

^^p{iuJn, ky) = £i°j(«^n, ka:, ky) + I]£i°J(^^n, k^, ky)^^{iuJn)^^p{iuJn, k^, ky) , (9) 
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with the summation over all planes, 7. The non-interacting planar Green's function, is similarly 
found by the Fourier transform 

G^liiuJn, kx, ky) = (^-^^ Q}°^ (^i' ^i' ^^n) exp [kx{xj - Xi) + ky{yj - yi)] . (10) 

The local Green's function, G(rj, rj, iujn), is required for calculating the self-consistent potentials, 
and the Green's function G{ri,Tj,iuJn) between two neighboring sites, rj = {xi,yi,Zi) and r^- = 
{xi, yi,Zi± a) is required for current calculations. These are given by the simple planar momentum 
integrals: 

G{ri,rj,iun) = {-] / / G JiUn,kx,ky)dkxdky (11) 

\aj J-TT/aJ-n/a ^''^ 



where again, a = Zi/a and j3 = zj/a, and the phase factors in the integral have canceled as Xi = xj 
and Hi = yj. Hence our goal is to find the interacting Green's functions, Q^pii^m k^, ky). 

We make a huge improvement (by one to two orders of magnitude) in the computational efficiency 
by transforming the two-dimensional planar momentum integral into a single integral over in-plane 
kinetic energy. In the case of nearest-neighbor hopping on a square lattice, the kinetic energy within 
the a-th plane is given by 

e^y = -2ta [cos(fc^a) + cos{kya)] = je^ (12) 

where ta is the hopping integral between two nearest-neighbor sites within the a-th plane and 
a is the lattice spacing. The effect of the in-plane kinetic energy is equivalent to an increase in 
the on-site energy, ef'^ h-* ef^ + £:*f which can vary between the different planes. The planar 
Green's functions only depend on the planar momentum via the normalized kinetic energy, = 
—2t[cos{kxa) + cos^kyo)], such that G.^f^{i^n, kx, ky) = Q_^^{iujn,e^y). Hence, by using the two- 
dimensional density of states, p'^^{e), for a square lattice, the momentum integral is transformed 
into 

/ ^ \ 2 p7v/a rn/a peso 

(tt) / G^p{iuJn,kx,ky)dkxdky= / G^p{iuJn,^)p^^(^)d^. (13) 

The specific hopping integral, ta, for each plane affects the on-site potential of a given plane in the 
continued fraction method [see Eq. (|I^) below], but does not contribute to the change of variables 
in the momentum integral. 

Once the system is converted to a one-dimensional model, with nearest-neighbor hopping, the 
Green's functions can be solved rapidly by a continued fraction expansion, without recourse to 
another fc-space integral for the 2;-direction. The continued fraction expansion is similar to the 
recursion methodS'i, modified to include superconductivityll, but with three important differences. 
First, the method is much faster, as there is no need to expand about a site to obtain a new basis 
— the system is already one-dimensional in form. Second, there is no inaccuracy in the termination 
process, as the hopping integrals are given exactly in the model. The third point is an alteration, 
because the sites of interest are not at the end of a chain, but in the middle. This leads to a different 
set of continued fractions that must be calculated compared with the standard recursion method. 
In test runs, our method proves to be 4 x 10^ times faster, and with machine precision accuracy, 
compared to a standard recursion method expansion which is terminated (due to memory limits) at 
an accuracy of one part in 10^! An alternate approach is to solve Dyson's matrix equation directly 
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for a finite system, where the infinitely extended bulk boundaries can be mimicked by appropriate 
choice of potentials for the end planes. We have carried out such an approach as a comparison, but 
find it to be much slower (by a factor of 4000 for 60 planes than our method, and it grows like A^^ 
for large systems with planes) so we only describe the continued-fraction method below. 

The equivalence of our method to the recursion method is that we calculate the Green's functions 
directly from a continued-fraction representation of the inverse of the Hamiltonian matrix in real 
space. That is, we find 

-1 







"■= =a—2 =a—l 



\ 







6^ , iuJnl — a , b 

=a—l =a—l =a 













=a =n+l 













b^ , 1 ioOn^ — Ci & , o 

=a+l "■= =a+l =a+2 



a[3 



where the matrices \ql^ are the total in-plane energies for a particular plane, given by 



and I^Q,} couple the a — 1th and ath planes. 



b 



'ta—l,a 

t 



(14) 



(15) 



(16) 



The local planar Green's function, G^^{iuJn,s^, is readily evaluated as a combination of continued 
fractions (as in the renormalized perturbation expansion^). We define the right-directed, R^{iuJn), 
and left-directed [i^n) continued fractions from a plane, a, recursively as 



iuOril — a — b , ,R , , (icu. 



iujnl — 



+1 



;i7) 

:i8) 



The recursive calculation continues to infinity, but once it has been extended to planes in the 
uniform bulk medium, where a < 1 — N/2 or a > N/2, the coefficients at each level become 
constant. The effect of a constant phase gradient in is equivalent to a constant phase factor in the 
hopping integral, ta,a+i, that does not change between planes in the bulk. Hence, by equating all 
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RjyiuJn,e''y) as R^{ioJn,e''y) for a > N/2 and L^{iuJn,e''y) as L^iii^n,^''^) for a < 1 - N/2 in the 
bulk limit, an exact terminator function can be calculated as the solution of a complex quadratic 
matrix equation: 



a - iunl] b^~^R {ioJn, e'^y) + b = 



a — iuj„ 1 
=00 "= 



b'^L (iujn,e'^y) + b^ = 0. 



(19) 
(20) 



Note that the same terminator function is used for all sites in the intermediate layers, and the 
functions R^ and calculated for one site are used in the calculation for the next site — so the 
number of computations required to find solutions for all sites can be less than 0{N). 

There are two ways to solve this matrix quadratic equation. When we perform calculations on 
the real axis, without including a supercurrent, the matrix equation becomes analytically tractable 
to solve. On the imaginary axis, we find that it is numerically faster to simply find an iterative 
solution to these quadratic equations within the bulk medium. In most cases, accuracies of one part 
in 10^° can be achieved in ten iterations or less. 

The continued fractions form the local planar Green's functions, according to 



G (iujn, e'^y) = \iiOnl - a (e'^y) - b^ , (icUn, e'^y^b - b R~l, (zu;„, e'^yW ] ^ 



(21) 



which, using Eqs.(|T^-[T^), can be simplified to 



Q^^iiuJn, e'^y) = \Rr,+ kry- ^^nl + 



(22) 



The Green's functions connecting neighboring planes, a and ail, which are required to calculate 
the current flow, are given in two equivalent forms 



lUJr. 



e-y) = -a^S'^n, e-y)k^^^R-l^{tuj^, e-y) = -L^\iuJn, £^^)^+i£,+i^„+i(^^n, e-y) (23) 

The local planar Green's functions enable us to calculate the self-energy and electron density, ?7,j, 
on each site in a given plane, a = Zi/a, with the latter given by 



/oo 
-oo 



(25) 



The current, Ia,a+i, which flows along each link between two neighboring planes, a and a + 1, in 
the 2;-direction is given by: 
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2,tc r°° 

= ^kBTj2 p^''{e-y)lm[Go.,o.+l{^uJn,e-y)]d^. (26) 

tl , 1 ^ —CO 

UJn 

A stringent convergence check for self-consistency, when there is a phase difference between the 
bulk superconductors, is that the current flow is constant from one plane to the next. 

For the bulk boundary regions, the uniform variation of phase in the off-diagonal self energy, 0(rj) 
has the form, 0(rj) = (pQ exp [zq ■ rj], where the net superfluid momentum depends on q = (0, 0, Qz) 
through mvs = ^ sm{qza) and leads to the following solution of Dyson's equation: 



and 



F{ri,rj,iun)= d-'k y^j^^v _ ^^28) 

J {lUn - + /i - {lUn + ^k-q - /i + - \(f)o{lUJn)\ 

where the diagonal self energy, T,{iu!n), is independent of site index, r^, in the bulk. We will need 
only the local, and nearest-neighbor bulk Green's functions in our calculations. 



C. Local self energy calculations 

In this contribution, the superconducting region is modeled by the negative-U Hubbard model 
within the Hartree-Fock (static mean-field) approximation. In this case, the local self energy is 
found from the local Green's functions by: 

T,{vi, iun) = UTJ2 Givi, Fj, iun) (29) 

and 

0(ri, iUn) = -UT J2 F{ri, r^, iu^), (30) 

where the instantaneous electron-electron interaction energy, U, leads to a time-independent self 
energy. This procedure is identical to the conventional Bogoliubov-de Gennes approachil, which 
neglects retardation effects in the superconductor. As all sites within a plane are identical, the self 
energy need only be calculated once for each of the planes. 

For planes within the barrier which include impurities, the dynamical mean field approximation 
says that the local (site) Green's function, G(rj, r^, iUn) is related to a local host Green's function^, 
g{ri,iun), via 
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(31) 



The atomic Green's function, (7°* (r,-, iuj„), which will be equated to the local Green's function, 
G_{ri,Ti,iuJn), in the dynamical mean field approximation, then satisfies 



G^jrj, iuJn) = (1 - Pimp)Q_{ri, iUn) + Pimp Q_ ^iju ^^n) " U^^X 



(32) 



and the local self energy becomes 



(33) 



Starting from a local self energy, ^(rj,iu;„), and a local Green's function, Q{ri,ri,iuJn), 
Eqs.( PT|P^J53| ) can be employed to iteratively determine a new self energy ^(rj,-ju;„) when the 
plane is described by the Falicov-Kimball model. The method, which is solved for a fixed concen- 
tration of impurities, pimp, is equivalent to the coherent potential approximation. The algorithm is 
summarized in Fig. 0. 



III. PHASE VARIATION 

Standard theory of Josephson junctions^ predicts the phase variation of the current in the weak- 
coupling limit to be I{6) = 1^ sin(^) where 6 is the phase difference across the barrier, and Ic is 
the temperature-dependent critical current. Such a current variation arises from consideration of 
two superconductors with different phases being coupled by a single energy-independent transmis- 
sion coefficient, corresponding to the tunneling of Cooper pairs across a barrier. A more general 
consider at ion0 includes the Bogoliubov-de Gennes equations for a one-dimensional system of two 
superconductors coupled across a potential barrier, V. The strength of the barrier is measured by 
Z = mV/h'^kF which determines the transmission coefficient, r, according to r = 1/(1 + Z^). It 
is found that Andreev bound states carry the current, whose variation with phase difference across 
the barrier leads toiil 

, = -G,^ , u„h i^^^pm] . (34) 

2e v'l-rsW(9/2) I iksT ) 



This reproduces the weak-coupling resultli^ in the limit of low transparency, r ^ 1, with 

— - sin(6') tanh ( , ° ^ , 
2e ^ ' \2kBTJ ' 



I = -Gn'^ sin(^) tanh f-^) , (35) 
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and reproduces the formula for a point contact Josephson junction with perfect transparency^, 
r = 1: 

j^-G,^,...(g/2)tanh( ^°;?y ). (36) 

e \ 2kBT J 

Interestingly, the maximum current is at a phase difference which approaches 6* = vr at low temper- 
atures in the latter case, whereas traditional tunnel junctions have their maximum at = 7r/2. 

In the more general case of a three-dimensional geometryU, the transmissivity depends on the 
angle of approach of the electron. The full current is obtained by an angular integral, da over the 
Fermi surface. The complete calculation leads to: 



^ ' 2e ^ ^ ^ J-./2 ^i-a^ia)sm\9/2) [ Sfc^T ) 

(37) 

where 

i?jv = / cr^vfa) cos(Q;)(ia. (38) 

J-tt/2 

The transmission probability for a single electron, which is proportional to its contribution to the 
conductivity, aj^, (from a Landauer formula) depends on its angle of incidence, as well as the 
strength of scattering at the barrier. 

The linear response of current due to a small phase difference, /', is given by: 

for both the one-dimensional and three-dimensional cases above. So if the barrier is parameterized 
by a single scattering potential, then the linear response supercurrent only depends on the normal 
state conductance according to Eq. ( |5D| ) and the product of /' ■ -Rat is a constant, independent of 
the microscopic details of the barrier. More realistic treatments of the barrier region are expected 
to result in deviations from Eq. (^) and allow I' ■ Rn to vary with the properties of the barrier. 

Further differences arise between methods which treat the barrier as a single scattering potential 
and our self-consistent treatment of the order parameter and off-diagonal Green's functions within 
the barrier. Our results indicate that the effective scattering strength of the barrier, Z, changes with 
temperature and current flow. In general, the proximity effect, which enhances coupling between 
superconductors, is weakened as the current flow approaches the critical current. Hence at large 
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current flow, the effective barrier is increased compared to its value in the linear response regime at 
close to zero current flow. The effect of self-consistency on the current-phase relationship appears 
to be more marked for weak barriers, where the current flow becomes relatively large. 

We quantify the current response by making use of the linear response I' = {dl /dO)^-^^^^, 
as well as the maximum current flow, J^, through the barrier. In the weak-coupling limit Ic is 
exactly equal to /', and our results show that in general the two are within 20% of each other, and 
scale almost identically with external parameters. As /' requires much less computational time to 
calculate, than Ic, we report values of I' for many of our results. 

IV. BOGOLIUBOV-DE GENNES RESULTS 

To begin, we demonstrate how our model reproduces standard results, by using the Hartree-Fock 
approximation to calculate the self energy within the barrier region. As such, we are effectively 
self-consistently solving the full Bogoliubov-de Gennes equations for the system. In all the results 
that we present in this paper, the superconducting region is modeled with an attractive Hubbard 
interaction of U = —2 at half filling. The homogeneous bulk superconductor then has a critical 
temperature, Tc = 0.11 and a zero temperature order parameter, Aq = 0.198. 

A. Barrier thickness 

We begin by solving for systems with a small attractive interaction within the barrier, Ub = —0.5, 
over a range of barrier thicknesses. This attraction is small enough that the bulk superconducting 
transition temperature of the barrier material is always less than any temperature we consider. 
Fig. ^ demonstrates the proximity effect, with the decay of the anomalous average at the center of 
the barrier as its width is increased. Fig. ^ is a log plot of the linear-response current and critical 
current against barrier thickness. As expected, both the linear-response current, J', and the critical 
current, Ic, drop rapidly when the number of planes within the barrier region is increased from 5 
to 30. The nearly constant slopes indicate that the decays are close to exponential. The exponent 
is approximately twice the bulk correlation length, given hj 2C, ^ hvf / A ^ 2ta/A ^ 10a, where a 
is the lattice spacing. 

Fig. 1^ shows the normalized current as a function of phase difference across the barrier 
region, for two different barrier thicknesses. The phase variation is compared to the sim- 
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pie form I{9)/Ic = sin(6'), which is appropriate in the weak-couphng hmit, and I{6)/Ic = 
(/o//c) sin(6'/2) tanh(A cos(^/2)/2/cbT), the result for a point-contact junction, appropriate in the 
limit of high transmissivity. The curve for a thin junction falls outside these two limits, indicat- 
ing that the effects of self-consistency and finite junction width are important in determining the 
current-phase relation for a Josephson junction. In particular, the calculated difference between a 
weakly coupled and strongly coupled junction lies in the opposite direction to the analytic formula 
on such a normalized curve. The actual magnitude of the critical current, is of course, greatly 
enhanced for the thinner barrier, as shown in Figs. 

B. Barrier interaction strength 

Fig. ^ demonstrates the proximity effect as a decay in the anomalous average. Fa = -F(rj, r^, r = 
0"*"), within a barrier region of twenty planes. The proximity effect is seen to depend on the 
Hubbard interaction, Ub, for sites within the barrier. The corresponding critical currents, Ic and 
linear response currents, are shown in Fig. |^. Note that in the cases where Ub is negative, the 
barrier region is actually a superconductor in its normal state, above its transition temperature. In 
the case where Ub = 0, the order parameter, Aj, is exactly zero within the barrier. In the example 
where Ub is positive, the order parameter actually switches sign within the barrier region, even when 
there is no external phase variation. The results of Fig. indicate that such a switching in sign of 
the pairing potential, Aj has no marked effect on the transport properties, which depend on the 
continuously varying Green's functions of the system. Once again we find a systematic tracking of 
I' and Ic with the strength of the Coulomb interaction, Ub, in the barrier. 

C. Barrier hopping integral 

Our method allows us to consider different hopping integrals in different regions (within planes 
or between planes). This would be appropriate when the barrier region has a density of states 
at the Fermi surface, or a Fermi velocity that differs from that found in the normal state of the 
superconducting regions. When modeling such systems, the hopping integral between successive 
planes that have differing intraplanar hopping integrals, is taken as the geometric mean of the two 
planar values = V^a-Wi)- 

Fig. ^ indicates that the critical current approximately scales with the hopping integral in the 
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barrier region at low temperature. This effect dominates over any increased scattering at the 
interfaces due to Fermi velocity mismatch. In fact, when th is large, the critical current is enhanced 
by a factor oi tb/t over that found in the uniform bulk system {tb = 1). The proximity effect results 
in a minimum of the anomalous average. Fa = F(rj,rj,r = 0+) =< Cj|(0^)cjj(0) >, at the center 
of the barrier region. At low temperatures. Fig. ^ indicates that the system with a smaller hopping 
integral in the barrier, tb = 1/2, has a larger anomalous average than the uniform system (with 
tb = 1), due to its increase in density of states at the Fermi surface. As the temperature is increased, 
there is a crossover, with the system that has largest hopping integral in the barrier {tb = 2), having 
the largest anomalous average just below the critical temperature, Tc. Such a crossover is due to 
the differing natural energy scales, tb, in the barrier region. For the system with tb = 1/2, an 
actual temperature of T = 0.1 corresponds to a temperature of 0.2 in the natural energy units of 
the barrier, tb- When the temperatures are given in units of tb, as shown in Fig. TH, a series of 
approximately parallel curves is seen, ordered according to the differing densities of states at the 
Fermi surface. The anomalous average. Fa, plotted on the y-axis, does not require such scaling, as 
it is dimensionless. 



D. Local density of states 

It is interesting to observe the density of states, in particular the presence of states within the 
gap as shown in Fig. [ll](a). The two major peaks correspond to Andreev bound states, which carry 
current, when there is a phase difference across the barrier region. The bound states can be seen to 
split in two, so that two inner peaks move towards each other, while two outer peaks move to the gap 
edge as the current flow increases. A careful examination of the curve corresponding to zero current 
reveals a great deal of structure, due to states trapped within the barrier region. The exact energies 
of the bound states depend on the width of the barrier, and the states arise in part from normal 
reflection of quasiparticles with non-zero momenta parallel to the planesSSEl. The intraplanar 
momenta of the quasiparticles may lead to the states being unobservable as current peaks at the 
corresponding voltages in an I-V curve, where only quasiparticles travelling perpendicular to the 



planes are measured. In Fig. |T2|, we examine the bound states in more detail, for a system with ten 
barrier planes (A^^ = 10). It is worthwhile noting that as the energy of the states within the gap 
approaches the gap edge, so they extend further away from the barrier, into the superconducting 
region [which begins after the fifth plane from the center [a > 5)]. The figure also demonstrates 
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an alternating parity between states. The states with energy closest to zero have even parity, as 
demonstrated by a maximum amplitude at the central planes {a = 0, 1). The next higher energy 
bound states have a node at the central planes, so exhibit odd parity. 

We can define a current-carrying local density of states such that the total current, Ia,a+i, 
between two successive planes, a and a + 1, is given by: 

Ia,a+1 = J ia,a+liE)dE. (40) 

A plot of the function ia,a+i{E) between the two planes, a = and a + 1 = 1, at the center of 
the barrier is shown in Fig. It can be seen that the majority of the current is carried by the 
Andreev bound states, and that the two peaks that have separated from a single Andreev peak at 
zero phase difference, carry current in opposite directions. The states at positive energy carry equal 
and opposite current to the states at negative energy, but their occupation is much lower for low 
temperatures, T < Tc. 



E. Charge Depletion or Accumulation Regions 

By including an extra repulsive or attractive potential on the interfacial planes which connect 
the barrier to the superconductor, we are able to model in a very simple manner, some of the 
effects of a charge accumulation region or a Schottky Barrier. In our simple model there is a mirror 
symmetry for a half-filled band : a repulsive potential depletes the electron density in one layer, 
to reproduce some effects of a Schottky Barrier, while an attractive potential results in a charge 
accumulation region. The effect on the pairing potential and current density are the same, for 
equivalent potential strengths at half-filling. That is, if a repulsive interfacial potential, Vi, results 
in a reduction in electron density on a specific plane, then the attractive potential, —Vi, at the 
interface, causes an equal in magnitude increase of electron density at that plane, and results in a 
system with the same variation in order parameter, and equal current response. 

Fig. |l5(a) indicates the oscillations in the anomalous average resulting from a single layer barrier 
with Vo = ±2t and Vq = ±4t. Notice that the proximity effect is reduced by the additional scattering 
potential, so the order parameter exhibits variation more like a step-function, with increasing barrier 
strength. We can think of the interfacial scattering as modifying the junction from SNS to SIS 
characteristics. 

Fig. depicts the effect on the electron density of a positive interfacial potential, which is like 
a Schottky barrier, since the electronic charge is depleted at the interface. The charge density also 
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exhibits Friedel oscillations away from the barrier. The transformation rii 2 — rii maps results 
for positive potential barrier to negative potential, Vq — Vq, and will change the system to that 
of a charge accumulation region. 



V. CHARGE IMPURITY SCATTERING 

We model a barrier region with impurities by using the Falicov-Kimball model, as described in 
Section II, using Eq. (||) with U^^ < and using the self-consistency procedure of Eqs. (|3Tl - |33D , as 



shown in Fig. |[ We carry out the calculations with an impurity concentration, pimp that ranges from 
0.01 to 0.2. In the limit of pimp = there is no scattering in the barrier, and the results correspond 
to the Hubbard model with = 0. Impurities in the barrier region lead to an imaginary part 
for the frequency- dependent self energy and, since the lifetime of the quasiparticles at the Fermi 
surface becomes finite, to a non-Fermi liquid, characteristic of annealed disorder scattering. 

A. Current Responses 



The addition of a small number of impurities is seen in Fig. |T6| to decrease the anomalous average 
in the barrier. The accompanying decrease in both critical current and linear response current, 
shown in Fig. |l^by the curves labeled UpK = —2, is more severe. Addition of 10% impurities {pimp = 
0.10) leads to a reduction in both current responses to approximately 1/3 of their initial values, while 
the anomalous average remains at approximately 3/4 its original amplitude. Anderson's theorem, 
which states that non-magnetic impurities do not detract from the superconducting properties of 
a system only holds for a spatially homogeneous system. A current flow breaks the symmetry, and 
a Josephson junction is inhomogeneous in one dimension, so the effects we observe do not violate 
Anderson's theorem. General considerations of Green's functions in a homogeneous system show 
that increasing the imaginary part of the electronic self-energy (hence reducing the quasi-particle 
lifetime) leads to a reduction in supercurrent for a given phase gradient. Hence it is expected that 
impurities would have a more deleterious effect on current responses than other superconducting 
properties, which is in the spirit of Anderson's theorem. 

We look at the linear response current, I', as a function of temperature for a barrier with 



impurity scattering, in Fig. |T^. The critical current remains slightly above zero until the bulk 



critical temperature, Tc. As suggested by Eq. (|39D, we can extract an effective transmissivity, 



18 



T = GN-h/2e^ of the junction from I' using: 

where Ao(T) is the value of the order parameter on the last superconducting plane before the barrier 
region. Fig. |19] depicts the results, which demonstrate that a real barrier has lower transmissivity 
with increasing temperature. This can be understood from the variation of the self-consistent order 
parameter, which is reduced to zero within the barrier as temperature is increased. These results 
are impossible to predict using more conventional approaches, and only arise from a microscopic 
model that includes self-consistency. 

B. Resistance Calculations 

We find self-consistent solutions of a system in the normal state, with no current flow, by 
setting the order parameter to zero on all planes. The self energy of planes outside the barrier 
contains only a constant real part, as initially we carry out the calculations within the Hartree-Fock 
approximation. Given the set of local self energies, the Green's functions coupling any two planes 
are readily found, for any momentum parallel to the planes. We are interested in the longitudinal 
components in the z-direction (perpendicular to the uniform planes) of the conductivity matrix. We 
define the conductivity tensor for our effectively one-dimensional system, from the linear current 
response Ia,a+i across a link between two planes, a and a + 1, due to an electric field, across 
all links between planes (3 and (3 + 1: 

(^a,l3 = (42) 

We find the conductivity matrix with frequency component, z/, neglecting vertex corrections (which 
is valid for homogeneous systems in the large dimensional limit) to be: 

(7a A^) = - / P^^i£xy)de^y [ ^ flm [Ga,/3+i{uj,e.^y)] Im [G^,«+i(a; + u,e^y)] + 

V \ ri J J-00 J-00 ZTT 

Im [Ga+iA^^ ^xy)] Im \Gp+iA^ + z/, e:r.y)] - Ihi [G , e^y)] Im [G^+i,„+i(w + z/, e^y)] - 
Im [G'„+i,^+i(a;, e^y)] Im [G^,a(w + e^y)]] [f{u) - f{u + v)] (43) 

where f{uj) is the Fermi-Dirac distribution function and p^^i^Sxy) is the two-dimensional tight- 
binding density of states, used for the sum over momenta parallel to the planes. 
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We are interested in the zero-frequency response, which is found from the appropriate hmit of 
Eq. (H): 

cosh^ {uj/{2kBT)) ^ ^ 

When calculating the resistance of the junction, it is important to be aware that for an inhomoge- 
neous ID system the current flow must be uniform but the electric field is not. The relationship 

Ia,a+1 = -^0 = X! ^a,f3Ei3^i3+i (45) 

leads to 

Efs,^^^ = j:{a-') lo, (46) 

a 

by multiplying on the left by the inverse of the conductivity tensor. The voltage across the junction 
is the sum of the electric fields across each link, multiplied by the lattice spacing, a, so we obtain 
the resistance, 

fl« = ^ (47) 

u a,l3 

given by the sum of components of the inverse conductivity tensor. 

It is worthwhile pointing out, that where there is no imaginary part to the self energy, so there 
is no quasiparticle decay, the conductivity tensor consists of constant elements, aa,[3 = o"o- In such 
a region, the electric field required to produce a current flow approaches zero as the inverse of 
the system size, so the local conductivity becomes inflnite. (That is, the sum of elements in one 
row of the conductivity tensor increases with system size.) However, the voltage drop across the 
region, given by the product of the electric fleld and length of perfect lead, remains constant (equal 
to lod/cTo, so the resistance is non-zero@ (equal to a/ao) while the local resistivity vanishes with 
large system size. In our calculations, where we neglect any lifetime effects of electrons outside the 
barrier region (in the Hartree-Fock approximation), there is still a contribution to the resistance 
of the junction from the 'perfect' leads, but the value of the contribution does not depend on the 
length of the leads, so it can be thought of as a contact resistance. 

Fig. |2y indicates the variation of junction resistance with impurity concentration in the bar- 
rier. The resistance increases linearly with number of scattering sites for the small concentrations 
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calculated, with the slope increasing with the strength of scatterers. The intercept is at a finite 
resistance, which corresponds to the resistance of an infinitely long, perfectly conducting lead, with 
conductivity tensor given by aa,f3 = ctq ~ 1.25. The resistance calculated for junctions with impu- 
rity scattering within the barrier region, does not change when the number of perfectly conducting 
planes on either side of the barrier is increased from one to twenty-five. Thereafter, numerical in- 
stabilities in the matrix inversion process make the calculations unreliable, but we can be confident 
that the answer already arrived at is the appropriate one for the infinite system. 

The product Ic-Rn decreases with increasing concentration of impurities in the barrier for the 
examples shown. That is, the reduction in critical current is greater than the increase in resistance 
due to impurity scattering. A system which increases the Ic-Rn product is found by incorporating 
an extra (coherent) superconducting region within the barrierEl. Such an SNS'NS or SIS'IS 
structure has been examined, where an SNS structure with 20 barrier planes has its central 6 
barrier planes replaced with superconducting material, creating a barrier sandwich of 7 normal, 6 
superconducting then 7 more normal planes. Both the critical current and linear response current 
increase by a factor of greater than two, while the normal state resistance is only reduced by 15% 
from its value for the SNS junction. 



VI. CONCLUSIONS 



We have developed an effective method for calculating the equilibrium properties of Josephson 
junctions. We are able to examine the microscopic details of self-consistently solved systems through 
local and non-local Green's functions. The importance of self-consistency has previously been 
shown@0, so we go beyond standard 'potential barrier' models of junctions to include the effects 
of spatial correlation and local fluctuations on the self-consistent potentials. 

Our results are interpreted in terms of a linear-response current, due to a small phase difference 
across the junction, as well as the critical current, being the maximum current that a junction can 
sustain. We flnd that a self-consistent microscopic determination of the potentials in the system 
for different phase differences results in current-phase behavior not predictable by standard flxed 
potential, semiclassical approaches. We study the effects of a number of properties of the barrier 
material, including its electron-electron interaction potential, its hopping integral and its charge 
impurity concentration. We also study interfacial scattering potentials, which mimic Schottky 
barriers and charge accumulation regions, and include barriers with a range of thicknesses. As well 
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as quantifying the change in current response due to such modifications in the barrier, we also depict 
the alterations in the proximity effect for all cases, and charge variation in the case of an interfacial 
potential. In future work, we plan to model the charge redistribution at interfaces more realistically, 
by incorporating a non-local Coulomb potential, which will self-consistently determine an effective 
potential with the charge density for each plane. By such a model, we expect to discover if there 
is any significant charge redistribution as a junction passes through its superconducting transition, 
as suggested by Greene's grou p@. 

We have also plotted the local density of states, to observe the bound states which occur at 
energies less than the bulk gap, within the barrier region. These states can be both current-carrying 
Andreev states, or states arising from normal refiections. Their evolution with position and phase 
variation is seen within our model, as is their contribution to the total current fiow when a phase 
difference is applied. We plan in later work to show how the detailed structure of electronic states, 
apparent in these equilibrium results, will affect the I-V characteristics of a junction. 

We have carried out resistance calculations for junctions with impurities in the barrier, and found 
that in general the reduction in critical current outweighs the increase in resistance due to charge 
impurities. We suggest that a junction of the form SNSNS, where a thin layer of superconductor is 
placed within an normal metal barrier, can increase the critical current of a junction dramatically, 
without markedly reducing its resistance. We are in a position to study how more subtle effects 
involving electronic correlation close to the metal-insulator can affect the Ic-Rn product. We plan 
to progress by making contact with specific materials, such as a Nb-InAs-Nb junction. 

Our microscopic formulation will be particularly necessary when we proceed to study junc- 
tions created with high-temperature superconductor materials. The d-wave symmetry of the order 
parameter, with its directional dependence, results in behavior that is not addressed by models uti- 
lizing a single transmissivity function. While an appropriate non-local version of the Bogoliubov-de 
Gennes equations can provide some insight into the properties of d-wave junctionsilS, use of the 
DCAMi is necessary to tackle the problem in any realistic manner. 
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FIGURES 




BULK 



FIG. 1. Microscopic stacked planar geometry of a Josephson junction. The sandwich of = 10 planes; 
Nsc = 4 superconducting planes coupled to a bulk superconductor on the left and = 2 barrier planes on 
the right, followed by a further Ngc = 4 superconducting planes coupled to another bulk superconductor 
on the right. The system is allowed to have spatial inhomogeneity only within the N modeled planes, but 
the calculations are always for an infinite system. 
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FIG. 2. Diagram of the iteration procedure, where the dynamical mean field theory is used to calculate 
local self-energies from local Green's functions. In the case of the Hartree-Fock approximation for the self 



energy, the dynamical mean field theory step is trivial, [see Eqs. (29-3 
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FIG. 3. Plot of the decay of the anomalous average due to the proximity effect in the barrier region. 
Note that the order parameter, A is equal to 2Fii = 2F(rj,rj,r = 0"*") in the superconducting region, 
where U = —2, and is equal to Fii/2 in the barrier region, as Ub = —1/2. It is the anomalous Green's 
function. Fa, rather than the order parameter, Aj, that is continuous throughout the system. 
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FIG. 4. Exponential decay of the linear response current, and critical current, Ic with increasing 
barrier thickness. Note how both /' and Ic systematically track with each other, and agree to within 20% 
for cases considered here. 
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Phase Difference, 9 

FIG. 5. Current flow as a function of phase difference across tlie barrier for a tliin and a thick junction, 
compared to analytic results for a low transmissivity (dotted line) and a high transmissivity (chain-dashed 
line) junction. Notice how the results for a thin junction lie outside of the results for the two analytic 
limiting cases. This shows how including self-consistency and a microscopic model of the barrier width can 
affect the current-phase relation. 
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FIG. 6. The decay of the anomalous average within the barrier region as a function of the Hubbard 
interaction, Ub- Note how the anomalous average decays more rapidly as the Coulomb interaction in the 
barrier, Uh, increases in value (and becomes repulsive). 
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FIG. 7. Linear response current, I', and critical current, /c, as a function of the Coulomb interaction 
in the barrier, U},. Note how both critical currents fall with U},, and note that there is no apparent 
discontinuity through [/;, = 0, where the sign of A changes within the barrier. 
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FIG. 8. Linear current, versus temperature for barriers with a hopping integral, tf,, that can differ 
from the hopping integral, t, between planes of the superconductor. Note how /' scales with ti, in the 
low-temperature regime. 
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FIG. 9. The variation with temperature of the anomalous average, Fa = F(rj,rj,r = 0"^), for the 
plane at the center of the barrier region. The different curves are results for different values of the barrier 

hopping integral, t^. 
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FIG. 10. The variation with temperature of the anomalous average, Fa = F(rj,rj,r = 0"^), for the 
plane at the center of the barrier region. The different curves are results for different values of the barrier 
hopping integral, tb, with the temperature normalized by th- 
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FIG. 11. The local density of states at the center of a Josephson junction (the plane a = 0), with 
increasing current flow through the device (Nj, = 20). The solid line is for zero current, the dashed 
line is for small current, and the chain-dashed line is for large current. The density of states within the 
superconducting region (plane a = 25) at zero current flow, where Aq = 0.198 is shown in the dotted line 
as a comparison. Note that the small bump at E = O.SAq for the superconducting region arises from the 
self-consistency relation and the proximity of this superconducting plane to the barrier. 
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FIG. 12. The evolution of the local density of states within the gap, as a function of position in a 
Josephson junction. The energy axis is normalized by Aq = 0.198. Note there are even parity bound 
states, with a maximum at the center plane, and odd parity states with a node at that point. The 
self-consistent solution is for a barrier of width 10 planes (AT^ = 10). 
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FIG. 13. The local current-carrying density of states, ia,a+i{E) at the center of a Josephson junction 
(where a = 0) and in the superconducting region (where a = 25). In this junction Nj, = 20. 
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FIG. 14. The off-diagonal local Green's function, Fa, indicating Friedel oscillations and a reduced 
proximity effect as the interfacial scattering is increased. 
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FIG. 15. Friedel oscillations in tlie electron density due to an intcrfacial scattering potential. 
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FIG. 16. The reduction of the anomalous average, Fa = F(rj,rj,r = 0"*"), for planes near the center 
of the barrier, with increasing impurity concentration, pj, and U^^ = —2. 
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FIG. 17. The critical current, Ic, and linear response current, I', as a function of charge impurities. 
Note that Ic and /' decrease as impurities are added to the barrier region. The greater scattering potential 
of U^^ = —2 results in greater decrease in current responses, compared to a scattering potential of 
UpK = —1- If superconducting planes are inserted in the middle of the barrier with U^^ = —2 (SNSNS) 
the critical current is seen to increase dramatically. 
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FIG. 18. The linear response current, /' as a function of temperature for increasing impurity concen- 



trations in the barrier. 
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FIG. 19. The effective transmissivity as a function of temperature. Note how the transmissivity decays 
with increasing temperature, shown here for barriers containing impurity scatterers. 
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FIG. 20. The normal state resistance of a junction versus impurity concentration in the barrier. The 
scattering potential is greater, = —2 for the higher resistance curve compared to U^^ = — 1 for 

the lower curve. Insertion of extra planes of superconducting material in the center of the barrier with 
U^^ = —2 (SNSNS) reduces the resistance in this case. 
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FIG. 21. The product of normal resistance of a junction with critical current, Ic-, and linear response 
current, for a junction with impurity scattering in the barrier. Scattering strengths of = — 1 

and = —2 are plotted. Increased figures of merit occur for both clean systems and for the SNSNS 

structures. 
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